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Abstract 


An all-at-once reduced Hessian Successive Quadratic Programming (SQP) scheme 
has been shown to be efficient for solving aerodynamic design optimization problems 
with a moderate number of design variables [1]. This paper extends this scheme to 
allow solution refining. In particular, we introduce a reduced Hessian refining technique 
that is critical for making a smooth transition of the Hessian information from coarse 
grids to fine grids. Test results on a nozzle design using quasi-one-dimensional Euler 
equations show that through solution refining the efficiency and the robustness of the 
all-at-once reduced Hessian SQP scheme are significantly improved. 

Key words, design optimization, constrained optimization, reduced Hessian meth- 
ods, quasi-Newton methods, successive quadratic programming, solution refining 

Abbreviated title. Reduced Hessian SQP with solution refining 


1 Introduction 

An aerodynamic design optimization problem can often be posed as 

min I(X,u) 

s.t. F(X,u) = 0, { ’ 

where X € W n denotes the discretized flow variables, and u € 9? m denotes the design 
variables, which, for example, could be geometry parameters describing the shape of a 
profile; I : K n+m — ► 3? is a cost function, which may, for example, measure the deviation 
from a desired surface pressure distribution; F : 9? n+Tn — ► §? n is a discretized version of the 
governing equations of the flow field. It is often the case that I and F are nonlinear, and 
the number of flow variables is much larger than the number of design variables (n >> m). 

In [1], a reduced Hessian SQP scheme is introduced for solving (1.1). This approach 
treats X and u as independent variables and updates them simultaneously at each iteration. 
One interesting property is that the flow equations are not required to be satisfied until 
convergence. It is intended to alleviate the cost of repeatedly solving the nonlinear flow 
equations required by other methods [9, 6]. Test results show that this scheme has a 
potential to be very efficient. This paper shows that the efficiency of the reduced Hessian 
SQP scheme can be further improved through solution refining. 

Solution refining techniques have been successfully used for aerodynamic design opti- 
mization by many authors [11, 7]. The basic idea is to solve the design problem on a coarse 
grid first, then use the coarse grid solution as an initial guess for a finer grid, and repeat 
this process until a solution on a desired grid is reached. Solution refining can be used with 
reduced Hessian SQP schemes as well. However in addition to refining solutions of flow 
and design variables it also needs to refine the reduced Hessian matrix approximation. This 
refinement is critical since the amount of Hessian information retained greatly influences 
the speed of convergence. However, as we are aware of, this issue has not been studied for 
aerodynamic design optimization. It is tempting to use the reduced Hessian approximation 
from the coarse grid solution directly as the initial guess for the next grid. However, this 
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practice may result in a mismatch of Hessian information between the two grids, particu- 
larly when unorthogonal bases are used for the null space, due to the dramatic change of 
bases in the null space caused by the change of grid sizes. 

This paper introduces a technique that satisfactorily solves the information mismatch 
problem. The basic idea is to use the fact that the reduced Hessian matrix can be treated 
as an invariant when the basis for the null space is properly chosen. This invariant bridges 
the information transition of the Hessian information between two grids. 

2 A review of an all-at-once reduced Hessian SQP scheme 
for aerodynamic design optimization 

Reduced Hessian SQP is a special case of SQP, which is a mature and successful technique for 
solving nonlinear constrained optimization problems due to its relatively low computational 
cost and fast convergence. 

In the context of aerodynamic design optimization problem (1.1), at each iteration of 
SQP, a quadratic approximation to the Lagrangian function of (1.1) 

L(X, u, A) = I(X, u ) + A t F(X, u), 

is minimized subject to a linearization of the flow equations. This gives a subproblem 

min d63? n+m g k d + \d T Hkd 
subject to F k + Ajd = 0, 

where d = ^ ^ ^ , g k = ^ |J ^ at (X k , u k ), A k = ^ ^ at (X k ,u k ), F k = F(X k ,u k ), 

and Hk is the Hessian (or approximation) of the Lagrangian function. 

One way of solving (2. 1-2.2) is via separation of variables. Suppose we can compute two 
matrices Yjt £ 3ff( n + m ) Xn and Zk E ^( n + m ) Xm such that the matrix [V* Zk] is nonsingular 
and Z%A h = 0 ( Zk is a basis for the null space of Aj). Let d = Ykd\ + Zkd<i and plug this 
into (2. 1-2.2). If we assume Z^H^Zk is positive definite (which is reasonable because this 
matrix is indeed positive definite close to the solution due to the second order sufficient 
optimality condition), by standard techniques, d\ and d<x are given by 

dx = -{A T k Y k )~ l F k (2.3) 

d 2 = -(Z k H k Z k )~ 1 (Z k g k + Z[ H k Y k d\). 

To avoid the computation of H k , the “cross” term ZjH k Y k d\ is simply ignored and 
ZjH k Z k is approximated by an m x m matrix B k using a variable metric formula such as 
BFGS, giving, 

d 2 = —B k 1 Z k g k . (2.4) 

Methods based on (2.3-2.4) are referred to as reduced Hessian SQP methods. 


( 2 . 1 ) 

( 2 . 2 ) 
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To avoid undue computational cost of orthogonal basis, many authors [3, 4, 2, 10] choose 
to use 


Z* = ( ) = ('/). (2.5) 

where, in aerodynamic design optimization, S is the sensitivity matrix. 

As analyzed in [1], the popular choice of Y* = [I 0] r potentially has an undesirable 
effect of producing a larger cross term. Instead 


( 2 . 6 ) 

is chosen, which is likely to produce a smaller cross term. Clearly ZjY* = 0 holds. 

From (2.3), d\ satisfies (AjYk)di = — F*, which is equivalent to 



which in turn is equivalent to 



J Jt (/ + 55 T )d 1 = -F fc , (2.7) 

where J k = at (X k ,u k ) is the current Jacobian matrix of the flow equation. (2.7) 

can be solved by first solving J k y = -F k for y, and then solving (/ + SS T )di = y for d\. 
The solution of the former is a linear flow calculation. The solution of the latter can be 
obtained by the conjugate gradient method, which is guaranteed to converge within (m + 1) 
iterations due to Rank(SS^) — m (see Golub and Van Loan [5]). Another way of solving 
(I + SS r )di = y is by inverting (/ + SS T ) directly. It is easy to show that 

(/ + SS T )~ 1 = / - S(I + S T S)~ 1 S T . 

Note that (/ + S T S ) is only an m x m matrix and its factorization can be obtained at 
minimal cost. 

After d\ and are available, d is given by 


d — Ykd\ 4* Zkd2 

= ( s T ) dl + ( / ) d 2 - 

Before going to the next iteration, we update the solution with 



3 



and the reduced Hessian approximation with the BFGS formula 


Bk + 1 = Bk - 


Bk£ks[]h VkVk 

sJ.BkSk Vk s k' 


where yk and Sk are given by 

Vk = ~ % k{ZjZk ) l Zlgk] 

Sk = (Zk+i^k+i)^ Zl+^ad, 


when certain update criteria are satisfied. 

In standard schemes the Lagrange multiplier is asked to satisfy 


(Y?A k )\k = ~ Yk9k , 


( 2 . 8 ) 


(2.9) 


( 2 . 10 ) 


which is equivalent to 

(I + SS T )(J T \ k ) = -Y? 9k. (2.11) 

In our scheme, Ajt is not explicitly needed, instead only the value of (J T Xk), which can be 
obtained by solving (2.11), has to be calculated. 

To ensure convergence, a merit function is needed to monitor the progress towards the 
solution. The l\ merit function is chosen for its simplicity and low computational cost, 
which is defined as 


^(X,tt) = /(X,u) + /x ||F(X,tt)||,. 

Throughout this scheme, the transpose of the flow Jacobian J T is never explicitly needed, 
which is desirable for many aerodynamic calculations. 


3 Reduced Hessian SQP with solution refining 

For the efficiency and the robustness reasons, solution refining is widely used in aerodynamic 
calculations. The basic idea is to get through the transient mode on cheaper, easier to solve 
coarse grids instead of on expensive, harder to solve fine grids. The translation of flow 
and design solutions from coarse grid to fine grid can be done in a variety of ways. Basic 
techniques include interpolation and mapping. Interesting discussions of related issues are 
covered by many papers including [11] and [7]. We will not go into details of them. Instead 
we will concentrate on the refining of the reduced Hessian matrix approximation, which has 
not been previously studied. 

For a given number of design variables m, an interesting fact is that the size of the 
reduced Hessian matrix is the same (m by m) on each grid level. On grid level /, assume (1.1) 
is solved with X[, u[, Z[ and B[. Through interpolation and mapping, an initial solution of 
the flow variables X[ +l and design variables Uj +1 at grid level 1 + 1 can be approximated. At 
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X' +1 ,tt' +1 on the new grid, the null space basis can be calculated. The question is how 
to obtain a reasonable level / + 1 reduced Hessian approximation at u l + l . One 

choice is B l m . However, B[ is intended to approximate the reduced Hessian Z l J while 
B l + l is intended to approximate the reduced Hessian Z[+ lT H[+ l Z[+ l . Since the Z matrices 
are not normalized, there is a potential mismatch of information between the two reduced 
Hessian approximations. A better choice is to treat the reduced Hessian in a normalized 
null space as invariant, which leads to 


(Z[ +1 (Z[ +lT Z[ +1 )~ l / 2 ) T H[+\z[ + \z[ +lT Z[ +1 )~ 1 ' 2 ) 


(zUzl'zlr^fHUzUz'Jzl)- 1 ' 2 ), 


which in turn implies 


Z' +lT fl' + 1 Z { +1 

* (z'+^z'+Y^zfzir^zftfiz: 

(Zi T Zi)- 1 / 2 (Z' +lT Z' +1 ) 1 / 2 . 


Hence it is reasonable to choose 


= (Z{ +lT Z{ + 1 ) 1 / 2 (Zi T Zi)- 1 / 2 Bi(Zi T Zi)- 1 / 2 (Z| +lT Z' +1 ) 1 / 2 . 


Since in our case a Z matrix always has full rank, Z T Z is always positive definite, 
and from its eigen-decomposition (Z T Z) -1 / 2 , as well as (Z T Z) 1 / 2 , can be obtained. Since 
the number of design variables m is much less than the number of flow variables n, the 
computation of the eigen-decomposition of Z T Z, which is of size m x m, is not significant 
and is affordable for many applications. In addition, since the decomposition is only carried 
out once at each grid level, the extra cost introduced by the Hessian refining is likely to be 
minimal. 


4 Quasi-one-dimensional Euler equations 

For our purposes here, we have chosen one popular form of central finite differences with 
nonlinear artificial dissipation, 

The quasi- ID Elder equations are 


F(Q) = d*E(Q) + H(Q) + D(Q) = 0 0.0 < x < 1.0 (4.1) 

where 


f p ■ 


pu 


0 

pu 

, E = a(x) 

pu 2 + p 

, H = 

-pd x a(x) 

€ 


,u(e + p). 


0 
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with p (density), u (velocity), e (energy), p = (7 — l)(e — Q.bpu 2 ) (pressure), 7 = 1.4 (ratio 
of specific heats), and a(x) = (1. - 4.(1 — a t )x( 1 — x)) (the nozzle area ratio), with a t — 0.8. 
For a given area ratio and shock location (here x = 0.7) an exact solution can be obtained 
from the method of characteristics. 

We choose one popular form of central finite differences to discretize these equations. 

d x q « 6 x qj = - — — j = 1, . . . ,Jmax 

Ax = 1.0/(Jmaar — 1), uj = u(jAx) (4*3) 

It is common practice and well known that artificial dissipation must be added to the 
discrete central difference approximations in the absence of any other dissipative mechanism, 
especially for transonic flows, see Pulliam[8]. 

For simplicity here, we use a constant coefficient dissipation of the form 

D 4 (Q) = V* (4.4) 


with 


~ = 9i+ 1 - Qj (4-5) 

with a typical value of = -j^j. 

Boundary condition at j = 1 and j = Jmax are defined in terms of physical conditions 
(taken from exact solution values) and will be treated as Dirchlet (fixed conditions) for now. 
The total system we shall solve is 


^(Q) = 


( 


tf J? E(Q) j - H(Q)j + Df(Q), j = 1, • • J N 
5 (Q)» = 0, i = 0, Jn 


(4.6) 


5 Test results for a nozzle design problem 

The nozzle design problem. We assume that a target velocity distribution u *, is given 
for each computational grid point. The design problem we are trying to solve is 

Find yi, i = 1, • • - ,m (spline coefficients describing a(x)J, such that 
\ Y2j=i X ( u j ~ u j) 2 ls minimized subject to (4-6) being satisfied. 

For our test examples, the breakpoints of the spline are evenly distributed in the interval 

[ 0 , 1 ]. 

Some implementation issues. At each grid level /, the design problem is solved using 
the reduced Hessian SQP to a specified tolerance e l (in our testing, t is uniformly set to 
l.E- 7), then grid is refined and the initial estimates to the design variables, flow variables, 
as well as the reduced Hessian on grid / + 1 are obtained from corresponding solutions on 
grid /. This procedure is repeated until the solution on the desired grid level is reached. 
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For our tests, a flat nozzle shape and a flat flow are used as starting data on grid 1. The 
starting reduced Hessian approximation is Bl = Zfzl 

Design results with grid size 65 and 7 design variables. In order to assess the 
efficiency of solution refining, we first choose to solve the design problem with grid size 65 
and 7 spline coefficients. This is quite a difficult problem because of over parameterization. 
We use the cost of solving the flow equations (without the presence of design variables) 
as a baseline to compare the cost of the reduced Hessian SQP with and without solution 
refining. 

Figure 1 shows 2 level flow solutions against the target solution. It also shows the 
convergence history of the flow residuals. The flow problem is first solved with grid size 33 
(upper subgraph). A significant number of iterations were taken due to difficulties in the 
transient regime (upper right). Then the solution is refined with grid size 65 (lower left). 
Due to good initial guess from the coarse grid solution, only 5 iterations were required to 
reach the tight tolerance (lower right). 

Figure 2 shows the design result of the reduced Hessian SQP without solution refining. 
Although the convergence of flow residuals and projected gradients appears to be superliner, 
it took a fair amount of iterations to get through the transient regime. The target, initial 
and final solutions are marked with solid, dashed and checked solid lines, respectively. 
The reductions in the projected gradient and the flow residual are marked with + and *, 
respectively. These notions are used consistently throughout the rest of this paper. 

Figures 3 and 4 show the test results of the reduced Hessian SQP with two levels of 
refining. First the design is solved on grid level 1 with grid size 33 (Figure 3). Then it is 
subsequently solved on grid level 2 with the desired grid size 65 starting from the grid level 
1 solution (Figure 4). Due to the good initial information supplied from coarse grids, it 
only took 10 (versus 35 without solution refining) expensive fine grid iteration to reach the 
tolerance. The overall saving over no solution refining in flop counts is about 33%. The 
final cost is only about 3 times as much as the cost of a single flow solution run. 

The test results are summarized in Table 1. 


Num. of Dsgn. Vars. 

Grid Levels 

Total Flops 

Num. of Iters. 

Cost Ratio 

0 

2 

2709573 

30(l)+5(2) 

1.000 

7 

0 

13631377 

35 

5.03 

7 

2 

9166365 

26(1)+10(2) 

3.38 


Table 1: Efficiency of solution refining (grid size 65) 


Design results with grid size 513 and 15 design variables. To make the design 
problem more challenging, we next solve the problem with grid size 513 and 15 design 
variables. This time the reduced Hessian SQP without solution refining failed to solve the 
problem (stuck in the transient regime), while the scheme with solution refining solved the 
problem successfully. 

Five levels of grid refining were carried out. Solutions on these five grid levels are shown 
in Figures 5 to 9. The results are summarized in Table 2. As expected, a relatively large 
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number of iterations were performed on grid 1 for getting through the transient regime 
(Figure 5). As the grid is refined, the convergence is getting smoother (Figures 6 to 9). On 
the final grid it only took 12 iterations to converge to the tight tolerance. 


Grid Level 

i 

2 



5 

Flops Count 

17570765 

15228630 



75681269 

Ml M II If I'HM 

39 

18 

14 

12 

12 

Grid Size 

33 

65 

129 

257 



Table 2: Results summary for grid size 257 and 15 design variables with solution refining. 

Testing with a various number of design variables and various grid sizes. Finally, 
we give the full set of the test results. The reduced Hessian SQP scheme with solution 
refining is tested on the nozzle design problem using a various number of design variables, 
i.e., 0, 1,3,7, 15 and various grid sizes, i.e., 33,65, 129,257,513. We want to point it out that 
each grid is associated with three flow variables. Hence the number of total flow variables 
is 99,195,387,771,1539 (3 times each grid size), respectively. Table 3 summarizes the test 
results, which is visualized in Figure 10. 

The problem with grid size 33 is solved with no solution refining, while the problem 
with grid size 65 is solved with 2 levels of solution refining, the problem with grid size 129 
is solved with 3 levels of solution refining, and so on and so forth. 

For each grid size, the cost with no design variables is used as a basis to measure the 
costs with a various number design variables. 

It seems that the cost of solving the design problem is about the same as or only one 
order of magnitude more than the cost of calculating the nonlinear flow, which is very 
promising for aerodynamic design optimization. Figure 10 shows that on the one hand, as 
the grid size increases, the relative cost of the design calculation is increased linearly. On 
the other hand it shows that as the number of design variables increases, the cost of design 
calculation appears to increase faster than linearly. We feel that this is caused by over- 
parameterization that results in much hard problems as more and more design variables 
are introduced for the same design problem. We believe that the increase should be linear 
without excessive over-parameterization. 


6 Concluding remarks 

This paper shows that solution refining can be combined with the reduced Hessian SQP 
scheme for aerodynamic design optimization. A particular useful technique introduced in 
this paper is the reduced Hessian refining under a variable metric formula. Test results show 
that the efficiency and the robustness of the reduced Hessian SQP scheme are significantly 
improved through solution refining. We believe that the efficiency could be further improved 
for certain problems through refining in the design variable space by gradually increasing the 
number of variables, particularly for problems with a significant number of design variables. 
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This brings forth a more challenging research topic for the efficient and effective reduced 
Hessian refining. Finally, applications of the reduced Hessian SQP scheme to 2D and 3D 
problems are currently under investigation by the authors. 


References 

[1] D. Feng and T. H. Pulliam, An all-at-once reduced Hessian sqp scheme for aero- 
dynamic design optimization, Tech. Rep. 95.19, Research Institute for Advanced Com- 
puter Science, NASA Ames Research Center, Moffett Field, California, 1995. 

[2] R. Fletcher, Practical Methods of Optimization, Second Edition, John Wiley & Sons, 
second ed., 1989. 

[3] D. Gabay, Reduced quasi-Newton methods with feasibility improvement for nonlinearly 
constrained optimization. Mathematical Programming Studies, 16 (1982), pp. 18-44. 

[4] J. C. Gilbert, Maintaining the positive definiteness of the matrices in reduced Hessian 
methods for equality constrained optimization, Math. Programming, 50 (1991), pp. 1— 
28. 

[5] G. H. Golub and C. F. V. Loan, Matrix Computations, The John Hopkins Univer- 
sity Press, second ed., 1989. 

[6] A. Jameson, Aerodynamic design via control theory, Journal of Scientific Computing, 
3 (1988), pp. 233-260. 

[7] G. Kuruvila, S. Ta’asan, and M. D. Salas, Airfoil design and optimization by 
the one-shot method, AIAA Paper 95-0478, (1995). 

[8] T. Pulliam, Artificial dissipation models for the Euler equations, AIAA Paper 85- 
0438, (1985). 

[9] J. Reuther, S. E. Cliff, R. M. Hicks, and C. P. van Dam, Practical design 
optimization of wing/body configuration using the Euler equations. AIAA paper 92- 
2633, 1992. 

[10] Y. Xie and R. H. Byrd, Practical update criteria for reduced Hessian SQP, Part I: 
Global analysis, Tech. Rep. CU-CS-753-94, Computer Science Department, University 
of Colorado, Boulder, 1994. 

[11] D. P. Young, W. P. Huffman, M. B. Bieterman, R. G. Melvin, F. T. John- 
son, C. L. Hilmes, and A. R. Dusto, Issues in design optimization methodology, 
Tech. Rep. BCSTECH-94-007 REV. 1, Boeing Computer Services, Seattle, Washing- 
ton, 1994. 


9 



Flow Solution 






Figure 1: Flow solutions with 2 levels of solution refining. Grid size = 33, 65. 
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solutions without solution refining. 7 design variables and grid size = 65. 


10 











1.2 


Flow Solution 


Nozzle Solution 



Optimal ity Conditions 





10 20 
Iterations 


Figure 3: Design solutions on grid 1 with 7 design variables and grid size = 33. 
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Figure 4: Design solutions on grid 2 with 7 design variables and grid size = 65. 
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Figure 5: Design solutions on grid 1 with 15 design variables and grid size = 33 
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Figure 6: Design solutions on grid 2 with 15 design variables and grid size = 65. 
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Figure 7: Design solutions on grid 3 with 


Flow Solution 




Figure 8: Design solutions on grid 4 with 
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15 design variables and grid size =129. 


Nozzle Solution 
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15 design variables and grid size = 257. 
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Figure 9: Design solutions on grid 5 with 15 design variables and grid size = 513. 
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Figure 10: Efficiency of reduced Hessian SQP scheme. 
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Table 3: Efficiency of the 
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Flops 

Cost Ratio 

2036491 

1.00 

2709573 

1.00 

4040497 

1.00 

6687365 

1.00 

11965583 

1.00 

1897820 

0.94 

3227728 

1.20 

5559882 

1.38 

10195908 

1.53 

20611283 

1.73 

2181578 

1.08 

4054968 

1.50 

7315226 

1.82 

13771675 

2.06 

28251307 

2.37 

5237019 

2.58 

9166365 

3.38 

16794298 

4.16 

29052365 

4.35 

56129543 

4.70 

17570765 

8.64 

32799395 

12.11 

55508501 

13.74 

93807404 

14.03 

169488673 

14.17 


1-at-once scheme 
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